POP Implementation of Mesoscale eddy diffusivity based on mixing theory
Here are my notes on the development of a parameterization for eddy diffusivity. The scheme is based on the paper by Michael Bates, Ross Tulloch, John Marshall, and Raffaele Ferrari. I'm working with Gokhan Danabasoglu and Matt Long here at NCAR and John Marshall (co-author) at MIT.
The general approach of this effort will be to follow John and Gokhan's implementation notes, isolating each of the terms contributing to the eddy diffusivity given by equation 6 in Bates et al. and comparing the CESM implementation of those terms with Bates, reanalysis, or other observations.
Implementation Notes on Diffusivity with Steering Level Suppression (Gokhan and John Marshall)
12 November 2014
Starting with \[K=u_{rms}∗L_{mix}\]
the general form for diffusivity, \(K\), is given by equation (6) of Bates et al. (2014). Namely,
\[L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} \]
Here,
\(u_{rms}\) is the root-mean-square (rms) eddy velocity;
\(L_{mix}\) is the mixing length;
\(L_{eddy}\) is the eddy diameter (depth independent);
\(u_{mean}\) is the mean zonal velocity (resolved);
\(c\) is the zonal eddy phase speed (depth independent);
- \(\Gamma = 0.35\)
\(b1 \sim 4\). ###\(\quad\quad\quad\quad\quad\quad\) Eddy Length Scales: Rossby deformation radius = \(L_r = {c_r \over |f|}\)
Equatorial Rossby deformation radius \(= L_{req} = {\sqrt{c_r \over 2\beta}}\)
Rhines scale \(= L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}\)
\(c_r\) is the first baroclinic wave speed computed following equation (2.2) of Chelton et al. (1998) with \(m=1\);
\(f\) is the Coriolis parameter;
\(\beta\) is the latitudinal variation of the Coriolis parameter; and
\(\sigma_{vi}\) is the Eady growth rate given by
\(\sigma_{vi} = {f \over \sqrt{R_i}}\)
with \(R_i\) the vertically integrated (over the 100 – 2000 m depth range) Richardson number.
So, any one these length scales could be used as an eddy length scale. An alternative is
\(L_{eddy} = min (L_r, L_{req}, L_{Rh}).\)
\(\quad\quad\quad\quad\quad\quad\) Eddy Velocity:
\({u_{rms} = alpha*\sigma_{vi}*L_r}\)
where \(\sigma_{vi}\) is the Eady growth rate based on local Richardson number and \(\alpha\) is a scaling constant.
\(\quad\quad\quad\quad\quad\quad\) Zonal phase speed:
\(c = - \beta * L_r^2\)
References:
Bates, M., R. Tulloch, J. Marshall, and R. Ferrari, 2014: Rationalizing the spatial distribution of mesoscale eddy diffusivity in terms of mixing length theory. J. Phys. Oceanogr., 44, 1523-1540, doi: 10.1175/JPO-D-13-0130.1.
Chelton, D. B., R. A. deSzoeke, M. G. Schlax, K. E. Naggar, and N. Siwertz, 1998: Geographical variability of the first baroclinic Rossby radius of deformation. J. Phys. Oceanogr., 28, 433-460.
Tullock, R., J. Marshall, and K. S. Smith, 2009: Interpretation of the propagation of surface altimetric observations in terms of planetary waves and geostropic turbulence. J. Geophys. Res., 114, C02005, doi: 10.1029/2008JC005055.
Questions:
In the 2D implementation, \(u_{rms}\) and \(u_{mean}\) specifications: upper-ocean vertically or integrated or at \(z = 0\)?
A: \(U_{rms}\) is not depth dependent; both \(u_{rms}\) and \(u_{mean}\) are for surface only
In the 2D implementation, vertical profile will be specified by \(N2(z)\)?
A: Yes, \(N^2(z) \over N_{ref}(z)\) to be more precise
Local \(R_i\) use imbedded in sigma in \(u_{rms}\) calculation?
A: No, \(u_{rms}\) is depth independent
alpha = ?
A: trial and error
Cancellation of \(f\)’s in \(u_{rms}\) calculation?
Zonal phase speed equation correct? Both \(\beta\) and \(L_r\) will be positive, >producing \(c < 0\) always. This appears to be in contrast with Tullock et al. (2009).
A: What is plotted in Bates is (U-c)
\(\color{green}{\quad CESM \quad Implementation \quad of \quad Steering \quad Parameterization}\)
During the development of the parameterization several tuning modifications were made to the terms listed above to better match the observations. * New implementation of Eady Growth rate (sigma) that avoids undefined values at the equator. See alternate derivation at end of this page. >\(\sigma_{vi} = {f \over \sqrt{R_i}}\) is replace by \(\sigma_{vi} = \quad {{\require\cancel f m^2} \over \cancel f N}\)
Don't use Rhine Scale in calculation of \(L_{eddy}\) >\(L_{eddy} = min (L_r, L_{req}, \cancel{L_{Rh}}).\)
Eddy phase speed calculation modified to avoid undefined value at the equator and limited to better match Bates paper.
>\(\cancel{c = - \beta * {L_r^2}}\) is replace with \(c = max(- \beta * L_{eddy}^2,-20)\)\(u_{rms}\) is hard to capture in this implememtation due to the coarse model resolution and use of mean vs resolved states. We imposed a \(5 cm/sec\) minimum on \(U_{rms}\) to tune to surface obs of Bates and included a scaling constant \(\alpha\) to handle the bias introduced by using mean state instead of resolved state >\(u_{rms} = max(5.,alpha*\sigma_{vi}*L_r)\)
\(\alpha=4\) The scaling constant for \(U_{rms}\)
- The Eddy diffusivity calculated by our implementation was weak so be bumped up \(\Gamma\) from .35 to 1.75 (5x increase) > \(\Gamma=.35\) replated with \(\Gamma = 1.75\)
\(\color{green} {Parameterizing \quad the\quad Eddy\quad Length\quad Scale} \)
NOTES: >\[K=u_{rms}∗{\Gamma * \color{red}{L_{eddy}} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} \] > >The Eddy Length scale is parameterized as \(L_{eddy} = min (L_r, L_{req}, \require{cancel}\cancel{L_{Rh}}).\) > >Since the two length scales we are using in the parameterization of \(L_{eddy}\) depend on the Baroclinic wave speed \(c_r\) we will first check \(c_r\) against Chelton 1998.
\(\quad\quad\quad\quad\quad\quad\quad\) First Baroclinic Wave Speed \(c_r\) CESM
\(\quad\quad\quad\quad\quad\quad\) Chelton Baroclinic gravity wave phase speed
The barclinic wave speed looks reasonable so lets compare the Paramaterized Rossby Radius to Chelton
\(\quad\quad\quad\quad\quad\quad\) First Baroclinic Rossby Radius CESM \(\quad L_{eddy}= min (L_r, L_{req}, \require{cancel}\cancel{L_{Rh}}).\)
\(\color{green} {Parameterizing \quad Zonal\quad Eddy\quad Phase\quad Speed}\quad (\color{red}{c})\)
NOTES:
\(K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - \color{red}{c}|^2 /u_{rms}^2 (z=0)} \)
\(\require{cancel}\cancel{c = - \beta * {L_r^2}}\) \(L_r\) too high at equator
\(c = - \beta * L_{eddy}^2\)
\(\quad\quad\quad\quad\quad\) Eddy Phase speed CESM
Compared to Hughes, the phase speed is much too high near the equator.
The first tuning mod is to limit the phase speed to 20cm/s.
\(\quad \quad \quad \quad \)Eddy Phase speed CESM \((c = max(- \beta * L_{eddy}^2,-20) )\)
\(\quad\quad\quad\quad\quad\quad\) Hughes Phase Speed (cm/s) from Tulloch Marshall Smith '09
\(\color{green}{Zonal\quad Velocity \quad Term}\quad (\color{red}{u_{mean}})\)
NOTES:
\(K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |\color{red}{u_{mean}} - c|^2 /u_{rms}^2 (z=0))} \)
Here \(U_{mean}\) is just CESM surface velocity.
\(\quad\quad \quad \quad\quad \) Zonal Velocity CESM
\(\color{green}{(U-c)\quad Term}\)
NOTES: \(c = max(- \beta * L_r^2,-20)\) )
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\) (U-c) CESM
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\) (U-c) Bates et al
\(\color{green}{(U-c)^2 \quad Term}\)
NOTES: \(c = max(- \beta * L_{eddy}^2,-20)\)
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (U-c)^2\) CESM
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (U-c)^2\) Bates et al
\(\color{green}{Parameterizing \quad Eady \quad Growth \quad Rate}\quad \color{red}{\sigma}\)
NOTES: >The original derivation of \(\sigma_{vi}\) is undefined at the equator because of \(f\)
\(\sigma_{vi} = {f \over \sqrt{R_i}}\)
The new derivation of \(\sigma_{vi}\):
$ R_i = {f^2 N^2 _}= $
\(\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}\)
\(\quad\quad\quad\quad\quad \)Eady Growth Rate \((\sigma) \quad\) CESM
\(\quad\quad\quad\quad\quad \)Eady Growth Rate \((\sigma) \quad\) Bates et al
\(\color{green}{Parameterizing \quad RMS \quad eddy \quad velocity}\quad \color{red}{u_{rms}}\)
NOTES: >\(K=\color{red}{u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} \) > >\(u_{rms} = alpha*{\sigma_{vi}}*L_{eddy}\) > >alpha (scaling constant) = 4 > >\(\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}\) > >\(u_{rms}\) is limited to 5 cm/s $ max(u_{rms},5.)$
\(\quad\quad\quad\quad\quad\quad\quad\) Eddy Velocity \((u_{rms})\quad\) CESM
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad u_{rms} \) Bates et al
\(\color{green}{u_{rms}^2 \quad Term} \quad\)
NOTES: >\(K={u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /\color{red}{u_{rms}^2} (z=0)} \) > >\(u_{rms}^2\) is a surface value
\(\quad\quad\quad\quad\quad\quad\quad\) Eddy Velocity Squared \((u_{rms}^2)\quad\) CESM
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad u_{rms}^2 \) Bates et al
\(\color{green}{{(U-c)}^2 \over {u_{rms}^2}} \quad \color{green}{Term}\)
NOTES: >\(K={u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * \color{red}{|u_{mean} - c|^2 /{u_{rms}^2} (z=0)}} \) > > >\(c = max(- \beta * L_{eddy}^2,-20)\) > >\(u_{rms} = max(u_{rms},5.)\) > >\(u_{rms}^2\) is a surface value
\(\quad\quad\quad\quad\quad\quad\) \({(U-c)}^2 \over {u_{rms}^2} \quad\) CESM
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad {|u-c|^2 \over u_{rms}^2} \) Bates et al
\(\color{green}{\quad Suppression \quad factor = {1 \over (1 + b1 * |\bar u - c|^2 /u_{rms (z=0)}^2 )}}\)
NOTES: > >\(c = max(- \beta * L_{eddy}^2,-20)\) > >\(u_{rms} = max(u_{rms},5.)\) > >\(u_{rms}^2\) is a surface value > >\(b1\) (scaling constant) = 4.
\(\quad\quad\quad\quad\quad\quad\quad\) Suppression Factor CESM
\(\quad\quad\quad\quad\quad\quad\quad\) Suppression Factor Bates
\(\quad\quad\quad\quad\quad\quad\quad\) Suppression Factor (Zonal x depth) CESM
\(\quad\quad\quad\quad\quad\quad\quad\) Suppression Factor (Zonal x depth) Bates
\(\color{green} {\quad Parameterizing \quad The \quad Mixing \quad Term \quad }\color{red}{L_{mix}}\)
NOTES: >\(\color{red}{L_{mix}} = \Gamma * L_{eddy} * Suppression\) > >\(Suppression= {1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2 )}\) > >\(\color{red}{L_{mix}} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}\) > >\(\Gamma = 1.75\) (Tuning mod: the original Gamma of .35 produced a Kappa with the correct structure but too weak.;
\(\quad\quad\quad\quad\quad\quad\quad\) LMIX CESM
\(\color{green}{\quad Parameterizing \quad Eddy \quad Diffusivity \quad} (\color{red}{K})\)
NOTES: >###\(\color{red}{K}=u_{rms}∗L_{mix}\)
\(\quad\quad\quad\quad\quad\quad\quad\) Eddy Diffusivity (K) CESM
\(\quad\quad\quad\quad\quad\quad\quad\) Eddy Diffusivity (K) Bates
\(\quad\quad\quad\quad\quad\quad\quad\) (Zonal x Depth) CESM
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\) K limited to (100 < K < 10000)
Here is the Bate's version
\(\quad\quad\quad\quad\quad\quad\quad\) Zonal average of the N2 normalized scaling CESM
The alternate derivation of \(\sigma_{vi}\):
$ R_i = {N^2 }$
$ N^2 = {{-g _0 }} $
After hydrostatic and geostrophic approximations
\(f \frac {\partial v}{\partial z} = {{-g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad f \frac {\partial u}{\partial z} = {{g \over \rho_0 }\frac {\partial \rho}{\partial y}} \)
so
\(\frac {\partial v}{\partial z} = {{-1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad \frac {\partial u}{\partial z} = {{1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial y}} \)
\(\therefore\)
$ R_i = {f^2 N^2 _}= $
\(\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}\)
\(RX_1 = RX_{east} = \Delta\rho_x = \rho_{i+1,j} - \rho_{i,j}\)
\(RY_1 = RY_{north} = \Delta\rho_y = \rho_{i,j+1} - \rho_{i,j}\)
\(RZ_1 = RZ_{k+1} = \Delta\rho_z = \rho_{k} - \rho_{k+1}\)
\(\displaystyle{1 \over L_{R_i}} \displaystyle\int_{2000m}^{100m} \left\lbrace { {-g\over\rho_0}{\frac {\partial \rho} {\partial z}} \over { {g^2 \over \rho_0^2 } \left[( \frac{\partial \rho}{\partial y})^2 + ( \frac{\partial \rho}{\partial z})^2\right] } } \right\rbrace dz\)
Note: missing \(f^2\) which will be cancelled when forming \(\sigma_{vi}\)
\(\quad\quad\) so \(\cdots\) this is not \(R_i\)
Implementation notes
Numerator : Top \(= -grav * RZ_{SAVE}(\cdots k+1) * dzwr(k)\)
Denominator :
$$
$$
\(work3 = {\left( TOP \over (grav^2*(work1+work2))\right)}*dzw(k)\)
Notes:
1)Need to be careful at top and bottom of ocean
2)Accurate dzw(k) for each (i,j) to form $L_{R_i}$
3) When constructing $sigma$ itself, use $RZ_{SAVE}$ with a minimum N value
4) use eps2